home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / data_smooth / smooth3.c < prev   
Encoding:
C/C++ Source or Header  |  1995-02-06  |  2.6 KB  |  85 lines
  1. /* Program 3. Smoothing and interpolation with any difference equation. */
  2. /* Contribution to Graphic Gems IV */
  3.  
  4. /* Paul H. C. Eilers, DCMR Milieudienst Rijnmond, 's-Gravelandseweg 565,
  5.    3119 XT Schiedam, The Netherlands, E-Mail: paul@dcmr.nl */
  6.  
  7. #define MMAX 100  /* choose the right length for your application */
  8.  
  9. typedef float vec[MMAX + 1];
  10. typedef float vecn[6];
  11.  
  12. void asmooth(vec w, vec y, vec    z, vecn a, float lambda, int m, int n)
  13. /* Smoothing and interpolation with any difference equation of order <=5.
  14.    Input:  weights (w), data (y): vector from 1 to m.
  15.    Input:  smoothing parameter (lambda), length (m).
  16.    Input:  coefficients (a) and order of difference equation (n).
  17.    Output: smoothed vector (z): vector from 1 to m. */
  18. {
  19.   static float b[MMAX + 1][6];
  20.   static int v[MMAX + 1];
  21.   int i, j, j1, j2, k, k1;
  22.   float s;
  23.   for (i = 1; i <= m + n; i++) {
  24.     v[i] = 1; if ((i <= n) || (i > m)) v[i] = 0;
  25.   };
  26.   /*  construct band matrix  */
  27.   for (i = 1; i <= m; i++) {
  28.     j2 = m - i; if (j2 > n) j2 = n;
  29.     for (j = 0; j <= j2; j++) {
  30.       s = 0.0; if (j == 0) s = w[i] / lambda;
  31.       for (k = j; k <= n; k++) s = s + v[i + k] * a[k] * a[k - j];
  32.       b[i][j] = s;
  33.     };
  34.   };
  35.   /*  compute cholesky-decomposition  */
  36.   for (i = 1; i <= m; i++) {
  37.     s = b[i][0];
  38.     j1 = i - n; if (j1 < 1) j1 = 1;
  39.     for (j = j1; j <= i - 1; j++) s = s - b[j][0] * b[j][i - j] * b[j][i - j];
  40.     b[i][0] = (s);
  41.     j2 = i + n; if (j2 > m) j2 = m;
  42.     for (j = i + 1; j <= j2; j++) {
  43.       s = b[i][j - i];
  44.       k1 = j - n; if (k1 < 1) k1 = 1;
  45.       for (k = k1; k <= i - 1; k++) s = s - b[k][0] * b[k][i - k] * b[k][j - k];
  46.       b[i][j - i] = s / b[i][0];
  47.     };
  48.   };
  49.   /*  solve triangular systems    */
  50.   for (i = 1; i <= m; i++) {
  51.     s = w[i] * y[i] / lambda;
  52.     j1 = i - n; if (j1 < 1) j1 = 1;
  53.     for (j = j1; j <= i - 1; j++) s = s - z[j] * b[j][i - j];
  54.     z[i] = s;
  55.   };
  56.   for (i = m; i >= 1; i--) {
  57.     s = z[i] / b[i][0];
  58.     j2 = i + n; if (j2 > m) j2 = m;
  59.     for (j = i + 1; j <= j2; j++) s = s - z[j] * b[i][j - i];
  60.     z[i] = s;
  61.   };
  62. };
  63.  
  64. void  pascalrow(vecn a, int n)
  65. /* Construct row n of Pascal's triangle in a */
  66. {
  67.   int i, j;
  68.   for (j = 0; j <= n; j++) a[j] = 0;
  69.   a[0] = 1;
  70.   for (j = 1; j <= n; j++) for (i = n; i >= 1; i--) a[i] = a[i] - a[i - 1];
  71. };
  72.  
  73. void gensmooth(vec w, vec y, vec z, float lambda, int m, int n)
  74. /* Smoothing and interpolation differences of order <=5.
  75.    Input:  weights (w), data (y): vector from 1 to m.
  76.    Input:  smoothing parameter (lambda), length (m).
  77.    Input:  order of differences (n).
  78.    Output: smoothed vector (z): vector from 1 to m. */
  79. {
  80.   vecn a;
  81.   int i;
  82.   pascalrow(a, n);
  83.   asmooth(w, y, z, a, lambda, m, n);
  84. };
  85.